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ABSTRACT: In this paper we describe a program (SOPHTY) implementing QED 
corrections to decays in the HERWIG++ event generator. In order to resum the 
dominant soft emissions to all orders, the program is based on the YFS formalism. 
In addition, universal large collinear logarithms are included and the approach can 
be systematically extended to incorporate exact, process specific, higher order cor- 
rections to decays. Due to the large number of possible decay modes the program 
is designed to operate, as far as possible, independently of the details of the decay 
matrix elements. 
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1. Introduction 

The use of Monte Carlo event generators has become an essential part of all ex- 
perimental analyses, both in interpreting data from existing experiments and in the 
design and planning of future experiments. Given the crucial role which Monte Carlo 
simulations play in experimental studies it is imperative that these simulations are 
as accurate as possible. 

While the existing Monte Carlo event generators have been highly successfully 
over the last twenty years, it was realised that a new generation of programs was 
necessary for the LHC. The reasons for this were twofold: firstly a number of new 
ideas to improve the accuracy of the simulations had been suggested e.g. [1-4]; sec- 
ondly the existing programs required major restructuring for long term maintenance 
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and to allow new theoretical developments to be incorporated. Given the changing 
nature of computing in high energy physics, the natural choice was to write these 
new programs in C++. A major effort is therefore underway, in preparation for the 
LHC, to produce completely new event generators [5], as well as new versions of 
established simulations [6-10] in C++. 

As part of the process of writing the new HERWIG++ event generator [8-10] we 
wish to improve many aspects of the simulation process. One area where improve- 
ments were needed was in the simulation of particle decays, both of the fundamental 
particles produced in the perturbative stages of the event and the decays of unsta- 
ble hadrons. Several major improvements have been made to the simulation of the 
decays: better modelling of the matrix elements in hadronic decays; and full spin 
correlations between the production and decay of particles [11]. Another problem 
with the particle decays in the FORTRAN version of the HERWIG program was the 
absence of QED radiation, which we will address in this paper. 

In existing Monte Carlo simulations the production of QED radiation in particle 
decays is normally simulated using an interface to the PHOTOS program [12-14]. 
This program is based on the collinear approximation for the radiation of photons 
together with corrections to reproduce the correct result in the soft limit [12,13]. Re- 
cently, it has been improved to include the full next-to-leading order QED corrections 
for certain decay processes [14]. 

Despite the success of PHOTOS it is based on the collinear approximation for 
photon radiation. The production of radiation in these decays is normally simulated 
in the rest frame of the decaying particle. The kinematics of many of the decays, 
particularly of the unstable hadrons, is such that the energy of the decay products is 
not significantly larger than their mass, in which case we do not expect the collinear 
limit to be a good approximation. However, there is always a soft enhancement 
for the emission of QED radiation. We therefore chose to base the simulation of 
QED radiation in HERWIG++ on the YFS [15] formalism for the resummation of 
soft logarithms. This formalism has the major advantage that the exact higher-order 
corrections can be systematically included, indeed the majority of the most accurate 
simulations including higher-order QED corrections are based on this approach [16- 
23]. 

Another significant improvement arises from the use of the C++ programming 
language and an object-oriented design for the program. The code framework gov- 
erning decays, in the HERWIG++ program, is arranged so that users can easily in- 
troduce the matrix element for a particular decay mode using the C++ inheritance 
mechanism. This framework also allows the inclusion of additional next-to-leading 
order matrix elements for both standard HERWIG++ and user defined decays. This 
makes it possible for the leading-order matrix elements and their higher-order correc- 
tions to be implemented in a systematic and consistent manor, rather than relying 
on one program to handle the leading-order decay and another for the higher-order 
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corrections. This will be of particular importance for the implementation of spin 
correlation effects. Moreover, it will be easier for users to introduce new decays, 
including higher-order corrections, as they will only require knowledge of, and make 
modifications to, one program rather than two. 

Note that, in full generality, it is not possible to consider radiative corrections 
to production and decay processes separately 1 , the minimal requirement for such 
a treatment to preserve gauge invariance is that the intermediate particle be on- 
shell. For most applications at hadron colliders we anticipate this to be a good 
approximation. At a technical level this amounts to neglecting off-shell effects in 
propagator numerators and including finite width effects in propagator denominators 
via an overall Breit-Wigner factor that multiplies the squared amplitude (the so- 
called narrow width approximation). This approximation is already in effect at the 
level of the tree amplitudes in the HERWIG++ generator, in which spin correlation 
effects are transmitted from the production stage to the decay stage according to 
the algorithm described in [1]. As the SOPHTY program is to dress tree-level events 
generated by the HERWIG++ simulation, a more subtle treatment including off-shell 
effects is beyond the scope of this work. 

The same approximation scheme is adopted in many other generators e.g. the 
PHOTOS and WINHAC [21] generators. Working in this approximation means that 
gauge invariance of the QED corrections is then guaranteed by considering only 
the universal leading log corrections to individualy decays in a cascade. Exact O (a) 
corrections, comprising of additional non-factorizable corrections are process-specific, 
and are therefore the subject of dedicated process-specific simulations; the SOPHTY 
paradigm is primarily one of universality and general applicability. Nevertheless, the 
dominant corrections are due to universal soft and collinear enhanced terms. 

In the next section we will present our master equation, based on the YFS 
formalism, for the generation of QED radiation, for the specific case of particle 
decays. High multiplicity i.e. greater than two body, particle decays are normally 
simulated as a series of sequential two body decays in HERWIG++. We therefore 
concentrate on the cases of the decay of a neutral particle to two charged particles and 
the decay of a charged particle to one charged and one neutral particle. In addition, 
we present algorithms for the event generation, using the master equation, for these 
two cases. The inclusion of higher-order corrections to the decays is then considered 
in section || followed by a discussion of the results of the simulation. Finally we 
present our conclusions and plans for further developments. 



1 The case of W boson production is a noteworthy example since it is both charged and unstable. 
This case has been studied in detail in the context of the charged Drell-Yan process [24-30] as well 
as single W production at LEP [31,32]. 
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2. Algorithms 



We begin by considering the n-body decay of a particle in the absence of photonic 
radiation. The partial width for such a decay is given by 

dTo = Wlj d$9 (2?r)4 S * ( P - S % ) P<# M « M *P ( 21 ) 
where 

d$ g = ff - (2.2) 

9 fi(2vr) 3 2gr 1 ' 

M is the mass of the decaying particle, p is its four-momentum and is the momenta 
of the zth decay product. We have also denoted the matrix element for the decay 
of a particle with helicity a by Ai a and p a p represents the spin density matrix. In 
order to simplify the expressions we have suppressed the dependence of the matrix 
element on the momenta of the decaying particle and the decay products. In practice 
we do not need to include the spin density matrix when calculating the total width, 
if we wish to average over the helicities of the decaying particles it is simply the 
identity matrix. However, inside the HERWIG++ simulation it allows us to include 
the correlation effects for the decay. 

In order to simulate the effects of additional QED radiation in the decay, we 



must generalise (|2.1|) to include the effects of arbitrary numbers of photons. In 
principle this extension is straightforward; one simply replaces the matrix element 
and augments the phase space. However, the matrix elements give rise to infrared 
divergences. The cancellation of these soft divergences must be made explicit before 
the Monte Carlo integration over the phase space can be performed. 

In the YFS formalism this cancellation of the divergences is manifest to all or- 
ders in perturbation theory. The cancellation relies on the fact that in the divergent, 
soft photon, limit both real and virtual corrections, to any process, take the form of 
universal kinematic factors multiplying the amplitude for that process without the 
additional radiation. In summing over all of the different soft photon contributions, 
these kinematic factors separately exponentiate, due to their universal nature. The 
resulting product of exponentials is the manifestly finite YFS form factor [15]. Resid- 
ual, non-factorizing, parts of the matrix elements, which cannot be exponentiated, 
are naturally infrared finite. 

Applying the YFS formalism to particle decays, by analogy to [17,18,33], we 



find that radiation modifies the n-body decay rate Q2.1]) to 

1 00 r 1 ™ 7 

r = E / d$ 777 eWC) II ^ PopMoMp C (2.3) 
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This is the master equation from which we intend to generate the QED radiation, 
where 

d$ = d$ p d$ fe (2tt) 4 5 4 (p - P - K) , (2.4) 

K = Y^a=i ki denotes the sum of photon momenta and P = Y^i=i Pi denotes the sum 
of the primary decay products momenta. The momentum of the ith decay product 
is pi (previously qi without the photon radiation), kj is the momentum of the jth 
photon, while d$ p and d$fc are the associated phase space measures: 

d% = Y\ f - ; (2.5) 
P fi(2vr) 3 2rf' 1 ' 

i=i ^ 

The symbol Q is used to denote the region of phase space inside which photons are 



< iv, with <v an energy 



soft and unresolvable, we choose this to be the region 
cut-off. We define 6 (k h Q) = for k { e Q and 6 (k h Q) = 1 for k t £ Q. This 
definition of f2 is not Lorentz invariant and in addition to specifying a value for iv 
we must specify the frame in which it has this value. 
The total dipole radiation function 

n 

Stotal {k)=Y,S{pi,Pj,k) (2.6) 

i<j 

is the sum of the individual dipole functions 

5 fe>R -,*) = ^ WA (_£i 5 -_S_) 2 , (2 . 7) 

where k is the four momentum of the photon, Zi is the charge of the ith particle 
in units of the positron charge and 6i = +1 (—1) if the momentum pi is outgoing 
(incoming). 

Likewise, the YFS form factor [15], Itotai {PiiPj, is given in terms of the form 
factor for pairs of charged particles 

n 

^totai (Pi, pj, n) = Y (p*>Pj> n ) • ( 2 - 8 ) 

i<j 

The YFS form factor for a pair of charged particles is given by 

Y ( Pi , Pj , n) = 2a (KeB ( Pi , Pj , Q) + B ( Pi , Pj )^ , (2.9) 

where Bij is the integral of the dipole radiation function over the soft photon phase 
space 

B( Pi , Pj ,n) = -L w (si f ^ (Pi- _ Pl_)\ (2.io) 

8tt 2 J n k \k-pi k-pjj 
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The virtual piece of the dipole 



does not depend on the cut-off, it is plainly Lorentz invariant. 

The form factors are given for the case of a neutral particle decaying to two 
charged final-state particles in Appendix [A.l| and a charged particle decaying to one 



charged and one neutral particle in Appendix |A.2| , in the rest frame of the decay 
products. For the case of a general moving frame only modifications to the Lorentz 
variant bremsstrahlung integrals are needed, in the final step where we have simplified 
with pi = —pj and Ei + Ej = M. The corresponding Bij functions are given in a 
general frame in [34] and [21] respectively. 

The factor C represents the total remainder of all of the matrix elements con- 
tributing to the total decay width for the particle, including any number of photons, 
when the leading soft divergent pieces are exponentiated and cancelled. The contents 
of C are referred to as the infrared residuals, they are infrared finite and are written 
as a power expansion in the electromagnetic coupling a. 

To order a there are three infrared residuals: the leading-order matrix element 
(O (a )), the finite remainder of the one loop correction to the leading-order process 
and the finite residual of the single photon emission matrix element. Using super- 
scripts to denote the order in a and subscripts to denote the number of emitted 
photons, we have 

C = 1 + o MM- [% (P ' {P ' ]) + ^ Pl t {P i'^ + ° ^) ' (2 ' 12) 
p aa >jvi a jvi a , y ^ 5 tota i(%) J 



where 



Po (P, {Pi}) = Pap {M a M l v * p + M\ a M% - 2aB total M a M*p) 
p\ (p, { Pi } , k) = PaP (i^M^Mfp - Stotai (kj) M a M%) 



(2.13) 



In (|2.13|) we use M.^ and My to denote the O (a) corrections to the leading-order 



matrix element (A4) from single real and single virtual photon corrections. The 
extension of the master formula to higher orders in the infrared residuals (j3) is 
straightforward, it is only limited by the usual technical difficulties associated with 
calculating Feynman diagrams that involve many loops or legs. 

In practice the leading-order matrix element is only strictly defined for the n- 
body phase space, not the phase space with additional photons. We therefore need a 
procedure which maps the momenta of the decay products after radiation, including 
the photons, on to the momentum configuration of the decay products prior to the 
generation of any radiation. Since the momenta of the primary decay products 
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are to be generated first before any QED radiation, according to the leading-order 
distribution, we can use these momenta to calculate the leading-order matrix element. 

In order to implement the theoretical framework of the master equation in an 
event generator, one expects that a number of different algorithms must be devised to 
cope with all possible decay multiplicities and charge configurations. In practice, due 
to the way in which particle decays are simulated in HERWIG++, most of the decays 
will either involve the decay of a charged particle to one charged and one neutral 
particle, or the decay of a neutral particle to two charged particles. Furthermore, 
we anticipate that more complicated decays will proceed via repeated applications 
and simple extensions of these two types of decay. We therefore concentrate on these 
cases. 



2.1 Final-Final Dipole 

The purpose of this algorithm is to dress the decays, generated by the core HERWIG+- H 
program, in which a neutral particle decays to two charged particles, with QED ra- 
diation. The input to our algorithm therefore consists of the momenta and quantum 
numbers of the parent particle and its children, distributed according to the leading- 
order differential decay rate fl2.1|) . 

In addition to infrared singularities, the dipole functions (|2.7| ) also exhibit mass 
singular behaviour associated with small-angle photon emission from the charged 
particles, in the massless limit. Therefore the angles of the radiated photons with 
respect to the dipole must remain fixed throughout the event generation process in 
order to produce an efficient, stable algorithm. 

In order to achieve this, we define the rest frame of the primary decay products 
and generate the photons in this frame. In this respect our approach is similar to that 
used in the /C/C event generator [20] for final-state radiation. Initially the photons 
are generated according to the dipole functions, which have a simple form in this 
frame. Implicit in the definition of the rest frame of the primary decay products (in 
this case the dipole) is the fact that the incoming three-momentum of the decaying 
particle must be equal to the total outgoing three-momentum of the photons. The 
three-momenta of the original decay products are then rescaled (reduced) to ensure 
energy and momentum conservation. 

A naive application of the method outlined above will lead to spurious results. 
It is important that we take into account the effects of the aforesaid choice of frame 
on the phase-space integration measure. To do this we employ the method of inte- 
gration over the Lorentz group, as described in [35], to transform the phase-space 
measure in (|2.3|). We start, by introducing the definition of the momentum of the 
decaying particle and the total momentum of the primary decay products in terms 
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of (^-functions, assuming the full phase space is to be integrated over i.e. 

J d$ = (2vr) 4 J d% d$ fc d 4 p ds d 4 P 2M5 3 (p) 5 (p 2 - M 2 ) (2.14) 

5 4 (p-P-K)5 4 [P-J^Pi ) S(P 2 ~s) . 



Secondly we insert the identity 



1 J V- 2 X lA X3 ( T-l P 



ix 7 s {-- l ) s [ L T») =l ' (2 ' 15) 

where L -1 is the boost from the frame in which X = ^X ,xj to the rest frame 
of X. This identity, in conjunction with those already present in ( |2.14j ), constrains 
the boost L to be the Lorentz transformation from the rest frame of the primary 
decay products (P) to the rest frame of the decaying particle (p) . We then change 
the integration variables, by applying the Lorentz boost L to all of the momenta 



involved, which is trivial as most of our expression ( |2.14|) is Lorentz invariant. This 

gives 

J d$ = (2tt) 4 J d$ p d$ fc d 4 p ds d 4 P d 4 X 5 3 (Lp) S (p 2 - M 2 ) (2.16) 

5 ^p-P-K)5 4 (p-±p)j5(P 2 -s) ^(v- 1 )* 
Integrating over the four momentum P, X and p we obtain 

J d$ = (2tt) 4 j d$ p d$ fc ds -^6 (M 2 - (P + K) 2 ) 5 4 (p - X>Y ( 2 - 17 ) 
The integral over s can then be performed giving, 

/d* = / d* r d* t -i ^(2w) 4 S l \P- V R ) . (2.18) 

As we first generate the momenta of the other decay products according to the 
leading-order matrix element we need to rewrite the integral in terms of the leading- 
order phase space. This is achieved by rescaling the momenta of the decay products 
before radiation to give the correct invariant mass for the decay products after the 
photon radiation. We define a momentum rescaling factor, u, such that the three 
momenta obey qi = uf>i. The momentum rescaling u is determined, by on-shell 
constraints, to be the solution of 



J2 \/ u2 \&\ 2 + m t ~ Vs = 0, (2.19) 
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where mf = pf = qf. 

The unitary algorithm techniques of [36] can be used to show that 

d* P s>( P -± P ) = / M s (q - ± „) n 



With this result we can rewrite our phase space measure as 



(2.20) 



r / n \ si i M -E?-i^ n /,o„3 



The decay width becomes 



(2.21) 



£ ^ / dr dd> fc e Wn)Jj5 total(fci)c ^ A ^JtII^ 



(2.22) 

with dr given by (|2Tl|). Equation ( |2.22j ) allows the construction of an algorithm in 
which the leading-order subprocess may be generated independently of and prior to 
QED radiation. 

Thus far we have treated the general case (an n body final state) but we will now 
specialise to the case of a neutral particle decaying to two charged particles. In this 
case the rescaling factor is u = \p\/\q\ where \p\ is the magnitude of the momentum 
of the decay products in their rest frame after the radiation and \q\ is the magnitude 
of the momentum of the decay products in their rest frame before the radiation. In 
this case the total width is 

£ / dr d$ fe J- e Wn) "Q £ {h) C^(l + ^\\ (2.23) 



where 



2M '' M' 

Up to now we have not made any approximations other than the truncation of 



the infrared residuals at (9(a). To simulate events using these results, ( |2.23| ) and 



(|2.22D , we need to make some approximations in order to obtain a distribution which 



is fast and efficient to generate by Monte Carlo methods. It is important to note 
that these simplifications are later exactly compensated by appropriate weighting and 
rejection of events. Naturally the first simplification we make is to omit the higher- 
order, infrared finite corrections represented by the factor C. We also neglect the 
factors associated with the rescaling of the leading-order phase space. Both factors 
tend to one in the limit of soft QED radiation and so neglecting them is reasonable, 
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given that the vast majority of photons produced will be soft. In addition to these 
two omissions we also approximate the momenta p\ and p 2 by the values they would 
have in the absence of any QED radiation, q\ and q 2 , this approximation is justified 
on the same grounds. These simplifications give the following crude distribution 

Tcrude = E / dr ° d$fe —f Y{qi ^ ] II $ (91, 12, h) . (2.24) 
n 7 =0^ 7 ' i=l 

Since we are working in the dipole rest frame, q\ = — q 2 , the kinematics simplify 
to the extent that the dipole function is analytically integrable. Moreover it means 
the only dependence of the QED part of ( |2.24j ) on q\ and q 2 is through their masses. 



Consequently we have the desired factorization that ( p.24| ) is really a product of two 
separate integrals, one for the leading-order decay and one for the QED radiation. 
The distributions may therefore be generated independently. Defining 

/d 3 k- 

yefefl) S{p 1 ,p 2 ,k i ), (2.25) 

we obtain 

X 

r CTUd e = r V — n n ~<e~ n = r . (2.26) 

/ — ' nJ. 

n~,=0 7 



In deriving (|2.26|) we have also made the approximation that the YFS form factor is 



Y m —n. According to r cru d e the photon multiplicity follows a Poisson distribution 
with average n. In practice it is sometimes useful to neglect part of the dipole 



distribution as described in Appendix B.l 



Once such a decay has been generated in the main HERWIG++ code it may 
dressed with QED radiation using the following algorithm: 

1. The number of photons is generated according to a Poisson distribution with 
average n. 



2. The momenta of the photons is then generated as described in Appendix |B.l 
This gives the crude distribution. 

3. The exact distribution (the master equation for Y) is obtained using rejection 
techniques. The weight for rejection is given by 

W = W dipole x Wyfs x W Jacobian x W higner , (2.27) 

where 



W — TT™ 7 g (Pl'P2» fc i) 

''Vdipole lli=l S{ qi ,q 2 ,ki)' 

w YFS = ^^1, 

Wjacobian = "^7^ + ^o) 

^higher — C. 



(2.28) 
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In practice the denominator of the dipole weight Wdi po ie is modified if we use 
the modified dipole without the mass terms. We denote a cut-off on the energy 
of the photons in the rest frame of the decay products as We will consider 
the contribution from the exact higher-order corrections in more detail in the 
next section. 

4. There is one remaining complication. The photons are generated in the frame 
where the decaying particle enters with momentum equal to the total photon 
momentum. However, we wish to apply the energy cut-off in either the rest 
frame of the decaying particles, or even the laboratory frame. There are a 
number of methods which we could use to achieve this. The simplest would be 
to evaluate the YFS form-factor in the rest frame of the decaying particle or 
the laboratory frame and veto any events in which any of the photons are below 
the cut-off in the relevant frame. However, this procedure can be inefficient if 
the veto removes a large number of events. 

We instead choose to use the same procedure as [20]. In this approach we 
neglect any photons which are below the energy cut-off in the relevant frame 
and apply an additional weight 

Wremove = exp (-Y 12 (q 1 ,q 2 ,tt B ) + Y 12 (qi,q 2 ,U) + Y 12 (pi,p 2 ,^) - Y12 (pi,p 2 ,^Ib)) , 

(2.29) 

where fl denotes the cut-off on the photon energies in either the rest frame of 
the decaying particle or the laboratory frame. 

For consistency in defining the infrared region, in applying this weight we do 
not apply dipole weights for those photons whose energy is below the infrared 
cut-off. 2 



2.2 Initial-Final Dipole 

In this subsection we describe our algorithm for dressing decays, in which a charged 
particle decays to another charged particle and a neutral particle. As in the final- 
final dipole case, the inputs to the algorithm are the momenta and quantum numbers 
of the parent particle and its children, distributed according to the leading-order 
differential decay rate. 

The situation here is less complicated than for the final-final dipole case because 
we can use the neutral decay product to absorb the recoil due to the photonic ra- 
diation. This allows us to simulate the radiation in the rest frame of the decaying 
particle. As with the final-final dipole we must also rescale the three-momentum of 
the charged particle in order to have overall energy-momentum conservation. 



2 In practice if we neglect the mass terms when generating the crude distribution we also need 
to include the weight from (B.8) for the removed photons as part of the dipole weight in order to 
take this into account. 



- 11 - 



As in the previous subsection, we begin by manipulating the phase-space mea- 
sure in order to factorize off a part of the integrand which can be interpreted as 
corresponding to the leading-order process. Taking p\ to be the momentum of the 
charged particle in the final state and integrating over the momentum \p[\ and p 2 , in 
the rest frame of the decaying particle gives 



2M ^ J nJ. 



where 



d$ 



1 



\\S {p, Pl ,ki) Pa pM a M^C (2.30) 



4(2tt) 

This can be rewritten as 



dfi Pl d$ fc 



Wit 



Pi(\pi\+ K cos 6 PlK ) +p° 2 \p 



(2.31) 



E 



n y =0 



dr d$ A 



\pi\ M 



1 _ e Y(P,Pun)Y[s(p, Pu k l )C 

cosO I + \q1\p2 \pi\ n " y ' 



l Qi\Pi ( \P 



K 



=1 



where 



dr r 



2^ d ^- 



?1 



;Paf3M a M* . 



(2.32) 



(2.33) 



4(27r) 2 M^' la " p - 

As before, by not changing the angles of the photons with respect to the dipole (in 
this case the charged final-state particle) we have dfl pi = dQ qi . The generation 
of the leading-order process (dr ) may proceed prior to, and independently of, the 
details of QED radiation. That is to say that no changes need to be made to the 
existing decay program, the QED algorithm for initial-final dipoles is universal in this 
respect. Momentum conservation and on-shell conditions require that the momenta, 
after generation of the photons, are given by 



p 



Pi = f V p 2 <zf + ml, pqtj , 
P2=[M- K 
K= (K ,K 



p 2 cp + mi, —K — pqi 



(2.34) 



where the rescaling factor p is 



A" 



COS#i#- 



((P 



1 +P2T + mf - ml) + (M - K ) J X (( P i +p 2 f ,m\,mi] ~ ±m\K\ 



2 1^| [(pi+ P 2r+Kl 



(2.35) 



and K 2 , 



K 



sin 2 6\k- 



The crude distribution, is obtained from the exact distribution (|2.33|) by drop- 
ping the kinematic factor arising from integrating over the delta function and the 
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higher-order non-soft photon corrections in C. The momenta in the form factor and 
dipole functions are replaced by the values generated from the leading-order decay 
(<fi = — q2) giving the crude distribution 

r crU dc =J2 dT ° d($>k —/ 12{mim J] ~ s (?» ?i> fc • ( 2 - 36 ) 

The dependence of QED part of ( |2.36| ) on the momenta q and q\ is overstated here, 
in the rest frame of the decaying particle the kinematics are so simple that this part 
only depends on the masses q 2 and q\. Therefore ( |2.36D is really a product of two 
independent integrals. The simplified kinematics allow the integral over the photon 
momenta to be performed analytically giving 

n = J ^9(^,0) S(q,q 1 ,k i ). (2.37) 



We therefore obtain 



crude 

=0^ 



oo 1 

V — n n -e- n = T . (2.38) 



In obtaining this we have, as in the final-final dipole case, approximated Y ps — n. 
Once again we have reduced the width to a simple Poisson distribution for the photon 
multiplicity. 

The generation of the crude width proceeds in the same way as for the final- 
final dipole. First we generate n 7 according to the Poisson distribution and then 
the photon momenta are generated according to the dipole functions (see Appendix 
|B.1|) . The form of the rejection weights W is similar to those in section [271] equation 
( [2.27D with the following changes: 

W _ TT™ 7 S{q,pi,kj) 

yVd\po\e - Hi=l gt q>quk .y 

w YFS = 



w \nlM (2 " 39) 

Whigher = C. 



Unlike the the case of the final-final dipole, we do not need a photon removal step 
because the decay is generated in the rest frame of the decaying particle. 



3. Higher-Order Corrections 

As stated in section ||, the effects of soft photons (photons with energy below the 
cut-off ui) have been included to all orders through the YFS form factor. If one 
neglects the infrared residuals in C, the effect of the master formula and algorithms 
is, for a given multiplicity, to generate the QED radiation according to the dipole 
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radiation functions only. This amounts to approximating matrix elements for the 
decay p — > pi...p n + n 7 7 by a product of eikonal factors multiplied by the leading- 
order matrix element. Ideally, we wish to include the higher-order effects in C as far 
as possible. 

Thus far our algorithms only require a set of momenta and their associated 
charges. Unfortunately calculating C exactly to a given order in a requires knowl- 
edge of the matrix elements for the specific decay process to that order. The structure 
of HERWIG++ is designed so that if these corrections are known they can be imple- 
mented. However, for the majority of decays these corrections will not be available 
and in this case we need to include the remaining enhanced contributions, i.e. the 
single collinear logarithmic terms. Depending on the mass scales involved, one can 
obtain a good approximation to C by just including these leading mass singular terms. 

In the collinear limit, the squared matrix element for a process including a mass- 
less emission, factorizes into the leading-order squared matrix element multiplied by 
a splitting function. The splitting functions only depend on the spin of the parti- 
cles involved. Therefore if, in addition to the momenta and charges, we supply the 
program with the spins involved in the decay, we may include the leading non-soft, 
collinear logarithms in C for the real emission contributions. 

In addition to affording us a way to include higher-order hard emission contribu- 
tions in a universal way, this approach has two further advantages. Firstly, as we shall 



describe in more detail in section |372| , in this approach we can readily obtain a good 
approximation to the virtual corrections. Secondly, given the logarithms associated 
to the collinear emissions are universal, they are necessarily gauge invariant. 

3.1 Real Emission Corrections:/^ 

In the quasi-collinear limit 3 , defined in [37], the matrix element including the emission 
of an additional collinear photon factorizes as 

n (=2^2 

i=l Pi ' K 

where M. is the matrix element for the leading-order process, \M\\ is the spin- 
averaged matrix element with the inclusion of one additional photon, Zj is the charge 
of the emitting particle, pi is the momentum of the emitting particle and k is the 
momentum of the emitted photon. Pa is the Altarelli-Parisi splitting function for 
emission of a photon from particle i, its form only depends on the spin of the emitting 
particle. 

In [38] these splitting functions were used, together with the factorization of the 
matrix element in the soft limit, to construct so-called dipole splitting functions for 



3 The quasi-collinear limit is the generalisation of the usual collinear limit to the case where the 
emitting particle is massive. 
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massive particles. These terms have the correct behaviour in both the soft and quasi- 
collinear limits and smoothly interpolate between the two, i.e. they reproduce the 
massive splitting functions for (quasi-) collinear emissions and the soft photon, dipole, 
radiation functions for soft emissions. We choose to use dipole-like terms based on 
the expressions in [38], omitting some sub-leading terms which were included in [38] 
to allowed the functions to be analytically integrated over the phase space of the 
emitted photon. With the dipole subtraction terms we may write an approximation 
for the real emission matrix element 



\ M r\ 2 ~ - e2 2 z & z i e i ( V v + %) \ M »\ 



(3.2) 



where indices i and j refer to the two particles forming the electric dipole and we 
have applied the conservation of charge ^2™ =o 0iZi = 0. 

We adopt the convention that the first index on TV, refers to the particle in the 
dipole which is considered to be emitting the photon, while the second index refers 
to the so-called spectator particle. From here we may write down a leading collinear 
approximation for the infrared finite residual (3\ 



Pi 



a 



An 2 



Paf3 M a M* p Z iQi Z fii [Ai + %]] 



(3.3) 



i<j 



where the are the infrared subtracted counterparts of V 



i.i 



Pi ■ k 



tyi ■ Pj 



(3.4) 



.(Pi+Pj)-k (Pi'k)_ 

For the case that both the emitter and spectator are in the final state, the dipole 
terms are given by 4 



V = -!- 



1_ 

Pi.k 
1 

Pi.k 
1 

Pi.k 



(Pi+P 3 ).k 



+ 



Pi.k _ 



(pi+k).pj (pi+pj).k 



2( Pj .k)(pi.pj) , 
a Pl +k). Pj ) 2 ^ 



2pj .k 
(Pj+k).pi 



Pi.k 
+ 



nr 



& 
Pi.k 



spin 0, 
spin ~, 
spin 1. 



(3.5) 



{Pi+P 3 ).k 

For the case that the dipole is comprised of the decaying particle (which we shall 
denote by index j) and one of its children, the Vij functions for emissions from the 
children are taken to be the same as in ( |3.5| ). However, for the decaying particle, we 
assume that it is sufficiently massive for us to neglect collinear enhancements, giving 
the following dipole function 



V ■ = -i- 



i_ 

Pj.k 



m j 1 
Pj .k 



spin 0, ~, 1 



(3.6) 



_(Pi+Pj)-k 

4 It is not possible to construct a quasi-collinear limit for the spin-1 splitting function, for a 
massive vector particle, in which the massless limit can be taken. We therefore use the corresponding 
expression for the massless dipole splitting function, augmented by a mass term which produces 
the correct behaviour in the soft limit. 
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Since the effects of collinear radiation from the decaying particle are neglected, only 
soft emissions are taken into account and so this dipole function is independent of 
the particle's spin. Furthermore, the infrared subtracted dipole splitting function, 
with the parent particle as the emitter, is identically zero: Dji = 0. 

In the soft limit these expressions reproduce the expected eikonal result 



and in the collinear limit the Dy equals the quasi-collinear splitting functions given 
in [37]. 

3.2 Virtual Corrections: /3q 

At present we have only implemented virtual corrections for two special cases in the 
SOPHTY code, those of initial-final and final-final dipoles with (relativistic) fermions 
in the final state, as in W and Z boson decays. These corrections turn out to have 
a negligible effect on distributions, compared to those of the real corrections. This 
is seen to be the case even for W — > e + z/ e and Z — ► e + e~ decays where one expects 
such effects to be greatest. 

As with the real emission we try to work in a universal way, without referring to 
the details of the matrix elements, using the leading log approximation. 

For both the case of the final-final dipole and the initial-final dipole the relevant 
virtual processes are represented by the lowest-order diagram with a photon joining 
the dipole constituents. On dimensional grounds, the large, leading logarithms of 
QED will be logarithms of M 2 /m 2 . Also, if we regularize the infrared divergences 
by introducing a fictitious photon mass m 7 , logarithms of M 2 /m 2 and m 2 /m 2 are 
possible. 

The infrared divergences from virtual corrections, must cancel the infrared di- 
vergences arising from the soft region of the photon phase space in the real emission 
process [39]. Likewise, terms diverging as m 2 — > 0, so called mass / 'collinear di- 
vergences, must also cancel between the real emission corrections and their virtual 
counterparts, this is the KLN theorem [40,41]. Using the fact that the m 7 — > and 
m 2 — > divergent logarithms have to cancel in this way, we can construct the leading 
log approximation to the loop integrals. 

To obtain the leading soft and collinear contributions to the virtual terms we 
therefore return to the soft and collinear approximation that was used for the real 
emission matrix element (|3.2|) . We calculate the leading logarithms arising in the 
real emission contribution by integrating the full dipole function over the full phase 
space for the emission of the photon i.e. both the soft k < uj and hard k Q > uj 
regions as was done in [38]. Performing the relevant integrals and taking the small 




(3.7) 
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mass limit gives the contribution of the virtual terms for the different types of dipole 

dF| LL = a (2 fin (M-) - l) In (£) + In 2 (g) + I In (g)) dr (initial-final), 

dF| LL = a (2 (in (S) - l) In (£) + iln 2 (g) + I In ($)) dr (final-final). 

(3.8) 

These expressions agree, at the level of large logarithms, with those obtained by 
direct calculation in [42] and [17]. From here we see that the /3q functions we require 
are, for the initial-final dipole 

« = (3 - 9) 

and for the final-final dipole, 

$ = -ln(^W (3.10) 

For resonant Z — > e + e~ processes this number is around 6%, dropping to around 3% 
for resonant Z — ► processes. The extension to other cases is obvious, it simply 

requires the use of the scalar and vector splitting functions instead. 



4. Results 

In this section we discuss the results from the SOPHTY program as implemented 
in HERWIG++. In order to test the algorithm we will consider both leptonic Z 
and W boson decays, due to their phenomenological importance. In addition we 
will consider a number of important meson decays to demonstrate the application 
of the program to hadronic decays. We reiterate that our approach simulates the 
soft photon corrections in the leading log approximation which depends on nothing 
more than the momenta and charges of the primary decay products, and simulation 
of hard collinear photons merely requires additional spin information. This being 
the case these examples represent a general test of our methods. 

The key distribution produced by the simulation is the total photon energy 
spectrum (K ). This is shown for Z decays in figure |l] and for W decays in figure 
|2|. We have considered a large range of masses for the decay products, including a 
fictitious heavy lepton (r'), to check for numerical instabilities and other irregular 
behaviour. For each type of decay we show the results of our algorithm including only 
soft photon effects and also with the dipole approximation for hard radiation. In all 
cases the amount of radiation is seen to decrease smoothly as the mass of the decay 
products increases, this can be understood from simple phase-space considerations. 

One can also see that the inclusion of the dipole splitting functions (T>ij) leads 
to an enhancement of hard photon radiation. This enhancement is most prominent 
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(a) Z Boson soft only (b) 1 Boson soft+collinear 
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Figure 1: The total energy (Kq) of the photons radiated in Z boson decays to leptons: (a) 
shows the Kq spectrum for the case that no infrared residuals are considered (C = 1); (b) 
shows the effect of including the collinear approximation for the O (a) residual 



(a) W Boson soft only (b) W Boson soft+collinear 
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Figure 2: The total energy (Kq) of the photons radiated in W boson decays to leptons; 

(a) shows the Kq spectrum for the case that no infrared residuals are considered (C = 1); 

(b) shows the effect of including the collinear approximation for the O (a) residual @\ . 



for lighter decay products, for the heavier decay products the effect is negligible. 
Again, this is to be expected as the mass of the emitting particle is known to screen 
the collinear divergence, this can be seen by considering, for instance, the massive 
splitting functions in [37]. 

Including the hard collinear enhancements also reveals a kink in the total photon 
energy spectrum. This kink occurs at a kinematic endpoint, beyond it all events must 
contain at least two photons which recoil against each other, hence the histograms 
drop beyond this value. Looking in this two photon region we also see that the 
photon multiplicity increases as the mass of the primary decay products decreases. 

Changing the spin of the primary decay products does not affect the soft distri- 
butions in figures ||a and ||a, it does however, influence the other distributions where 
hard collinear photon effects are introduced. The program uses the other splitting 
functions in ( |3.5| ) to account for this, although in the case that the decaying particle 
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(a) WINHAC vs SOPHTY (b) WINHAC vs SOPHTY: Relative Difference 




O 10000 20000 30000 40000 50000 10000 20000 30000 



Figure 3: The total energy (Kq) of the photons radiated in W^ 1 — > e v e /v e decays. In 
figure (a) the red histogram was generated using the WINHAC [21] simulation, including 
exact O (a) real emission corrections to the W^ 1 decay, while the black line was generated 
using the SOPHTY module for QED radiation in HERWIG++. In figure (b) we show the 
difference between the spectra shown in (a) divided by the associated statistical error. The 
discrepancy in the region beyond 40 GeV is exclusively comprised of events with at least 
two non-soft photons, which neither program is designed to model well. 



is a scalar there will be no collinear enhancement since in this case = 0. 

In figure f| we compare the total photon energy spectrum for W — ► ev e decays 
as generated by our program and that of the WINHAC generator. The agreement is 
seen to be quite good except in the region beyond the kink at around 40 GeV. As 
mentioned earlier, this region is populated exclusively by events with at least two 
hard photons. Consequently neither simulation expects to model this accurately. A 
correct modelling of this region will require the calculation of the infrared residuals 
(C) to O (a 2 ). This extension may be implemented in future versions of the program. 
We note that WINHAC was been independently compared with another simulation 
of the charged Drell-Yan process, HORACE, in [43], where good numerical agreement 
between the different approaches to QED radiation was recorded. 

The total energy of the photons radiated in the decays of some neutral vector 
mesons is shown in figure [|. Here we see that the energy distribution shows a 
behaviour that qualitatively resembles that seen in the case of the Z boson. In 
this the decay products are pseudoscalar mesons, there is no effect from 

including the collinear approximation for the radiation. In addition figure |]c shows 
the radiation for an example, K*° — ► K w T , with unequal masses for the decay 
products. Here we see that the lighter decay product is responsible for the more 
energetic photons, as we increase its mass (for illustrative purposes) the distribution 
tends toward lower energies. 

The total energy of the photonic radiation in the decays of some charged vector 
mesons is shown in figure [5|. As expected, for these decays the distributions be- 
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(b) Decays 
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Figure 4: The total energy (Kq) of the photons radiated in the decays of neutral vector 
mesons to pseudoscalar mesons for a number of different decays: (a) p — > 7r + 7r~; (b) 
4> vr+vr" and K+K"; (c) K*° -> K ± vr =F ; (d) J/-0 -> tt+tt" and J/ty -> K + K~. In 
addition to using the real physical masses of the decay products we have included the effect 
of varying the masses of the decay products. 

have in a similar way to those of the W boson, since they involve the same type of 
dipole. As with the neutral vector meson decays there is no effect from including 
the collinear approximation for the photon radiation as the decay products are pseu- 
doscalar mesons (T>ij = 0) . The K*^ 1 decays show the effect of having unequal mass 
decay products. 

For the leptonic decays of the charmonium resonances and the T resonance, the 
total energy spectrum of the radiated photons is shown in figures | and [7| respec- 
tively. As for Z decays, the effects of varying the r mass are included to show the 
mass dependence of the results. In the charmonium decays to t + t~ pairs we see 
a suppression of QED radiation since these decay modes are near the production 
threshold; for the physical r mass, this decay mode is not accessible in J /if) decays, 
while for if) (2s) and if) (3770) it is just below the threshold. The T resonance is sig- 
nificantly more massive and therefore the associated photon energy spectrum more 
closely resembles that seen for the case of the Z boson. 

These results show that the approach can be successfully applied to both per- 
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Figure 5: The total energy (Kq) of the photons radiated in the decays of charged vector 
mesons to pseudoscalar mesons for the decays: (a) p ± ti^tf ; (b) K** K ± vr° and 
K*^ 1 — > K ^. In addition to using the real physical masses of the decay products we have 
included the effect of varying the masses of the decay products. 

turbative decays and non-perturbative hadronic decays. 



5. Conclusions 

In this paper we have presented a universal theoretical framework for calculating 
QED radiative corrections to particle decays based on the YFS formalism [15] and 
the methods of [17, 18,33]. The essence of this approach is a reorganization of the 
perturbation series to resum all soft divergent QED logarithms. This formalism led 
to the master formula presented in ( |2.3|) . 

The master formula forms the basis of the Monte Carlo event generator SOPHTY, 
which provides QED radiative corrections for decays inside the HERWIG++ gener- 
ator. The Monte Carlo simulation takes into account large soft photon logarithms 
to all orders. In addition, the leading collinear logarithms are included to O (a) 
by using the so-called dipole splitting functions and inferring the associated virtual 
corrections with the aid of the Bloch-Nordsieck and the KLN theorems. 

Algorithms for the two basic "building-block" cases of dipoles comprising either 
two final-state particles, or the initial-state particle and one of its decay products, are 
presented in section |[ Although integrals like that of the master equation can gen- 
erally be readily performed with conventional Monte Carlo methods to give weighted 
events, the manipulations required in order to produce unweighted events with good 
efficiency (i.e. an event generator) are non-trivial. 

In designing these algorithms we were constrained by the requirement that the 
QED radiation should be generated, as far as possible, independently of the details 
of the main HERWIG++ program, which should provide the initial distribution of 
decay products. This was made possible due to the form of the master equation, 
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(c) ^(3770) soft only 
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Figure 6: The total energy {Kq) of the photons radiated in leptonic decays of charmonium 
resonances is shown above for J/ip (a/d), tp(2s) (b/e) and -0(3770) (c/f). The distributions 
on the left, figures (a), (b) and (c), are obtained by truncating the infrared residuals at 
0(1), whereas in (d), (e) and (f), the dipole splitting functions are used to include the 
effects of hard collinear photons. In addition to the real charged leptons we have included 
the effect of varying the r mass to illustrate the mass dependence of the results. 

the universal nature of the radiative corrections involved and our simplified crude 
distribution from which we initially generate the photons. Care was taken to design 
the algorithms to keep event weights as close to one as possible and to avoid numerical 
instabilities. Key to realising these features are importance sampling techniques and, 
importantly, a careful choice of frame in which to generate the radiation. 
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(a) T soft only 
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(b) T soft+collinear 
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Figure 7: The total energy (Kq) of the photons radiated in leptonic decays of the bottomo- 
nium resonance T. The distribution on the left (a) is obtained by truncating the infrared 
residuals at 0(1), whereas in (b), the dipole splitting functions (Pjj) are used to include 
the effects of hard collinear emissions. 



Our algorithm was tested successfully for several different types of particle decay 
produced by the HERWIG++ generator. In section f| we have shown results for 
the important cases of W and Z decays to various lepton species. In all cases the 
distributions show a smooth and stable behaviour agreeing with our expectations. 
A preliminary comparison of the total photon energy spectrum from the WINHAC 
generator shows good agreement and provides a good check of our methods. The 
application of the program to both hadronic and leptonic meson decays was also 
illustrated in section |j. 

There are several possible extensions of this work, for instance, there are a num- 
ber phenomenologically important decays for which the full O (a) corrections are 
known e.g. W and Z decays. As the code is designed to readily allow these cor- 
rections to be included they will be implemented in the near future. In addition, 
there are a small number of cases inside the HERWIG++ where we have to simu- 
late multi (i.e. greater than two) body decays and the extension of the algorithm 
to these cases would be useful. This is the first use of the YFS approach within 
the HERWIG++ program and there are a number of other potential applications, 
for example initial-state radiation in lepton collisions, which may be pursued in the 
future. 

In conclusion, we have applied the YFS formalism to the simulation of QED 
radiation in particle decays. The simulation, SOPHTY, based on the results of this 
work can simulate QED radiation in a wide range of particle decays and will be 
available in the next version of the HERWIG++ program. 
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A. YFS Form Factors 

In this appendix we give the expressions we calculate for the YFS form factors. For 
both the final-final dipole form factor and the initial-final dipole form factor we use, 
Pi = \pi \ /Ei , to denote the velocity of particle i. Furthermore, in each case we have 
assumed that the momenta obey p = Pi + P2- 



A.l Final-Final Dipoles 

Below in (|A.1|) we have the YFS form factor for a pair of final state charged particles 



with momenta pi, P2 whose combined momentum is p. Expression (|A.1|) is valid in 
the frame p\ = —p2- 
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We have checked that for the case /3j 
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in agreement with the previously obtained results [19]. 



A. 2 Initial-Final Dipoles 

In equation (|A.3| ) we have the YFS form factor for a pair of charged particles of 
momentum p and pi, (p2 = p — Pi), evaluated in the rest frame of p. 

Y(p,p 1 ,n) = %z p z 1 Y(j>,p 1 ,n) 

Y (p, Pl , a) = In (£) + In (g) - i In (f±| 

- ^ln 2 

2ft 
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We use Z p to denote the charge on the particle with momentum p. For the special 
case that (p — p 2 ) 2 is zero, as in leptonic W ± decay with a massless neutrino, we use 
the a more specialized compact form, to avoid potential numerical problems: 

Y( P , Pl ,n) = % z p z x Y( P , Pl ,n) 

f(p, Pl) Q)=ln(^)+ln(^)-iln(i 



2ft ; V 2ft / \ 2ft / V 2/3 




We have checked, analytically, that in the limit fa — > 1 the virtual contributions 
(He B (p,pi)) to Y (p,pi,Q) inside (|A.3|) are equal to those in (|A.4j) . In both cases 
the real contributions are identical, they do not involve fa, naturally as these con- 
tributions should only involve the moving charge in the dipole. Finally, as a check 
we observe that, dropping terms smaller than O (p\/p 2 ) the form factor ( |A.4|) agrees 
exactly with the corresponding expression in [21]. 



B. Generation of the Dipole Distributions 

In this appendix we describe how to generate the photon momenta from the dipole 
radiation functions. 
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B.l Final-Final Dipoles 

Consider the integral of the dipole function in the rest frame of p\ + p 2 

[ ^ S( Pl ,p 2 ,k) = -ZL Zl Z 2 [ dc d0 dA; k - -^rV- (B.l) 

J k An 2 J \pi ■ k p 2 ■ k) 

We choose to define the photon momenta as being with respect to p±, with c = cos 6* 
i.e. 

Pl .k = f?ifc (l -Pic) , , B2 ^ 
p 2 .k = E 2 k (l+(3 2 c). [ ' ' 

Using this representation of the momenta the integral can be rewritten as 

f*W S(pi,p 2 ,k) = -ft Zl Z 2 Jdc d0 dlnA; + {1 %%% 2C) - ^) . 

(B.3) 

The photon momenta can be generated according to this distribution in the following 
way. 

1. First the magnitude of the photon's momentum is generated logarithmically 
between cu, the minimum photon energy cut-off, and the maximum possible 
photon energy, E max , i.e. 

E x n 

-'-'max 



k = coi^ , (B.4) 

where TZ is a random number uniformly distributed between and 1. As the 
photon momenta are generated in the rest frame of the dipole the maximum 
energy of the photon is 

£ „ = | (_*_ _ -±-) . (B.5) 
2 \mi +m 2 M J 

2. The azimuthal angle is randomly generated between and 2n. 

3. The generation of the polar angle 9 is more complicated. The polar angle is 
generated by only using the interference term i.e. neglecting mass terms. This 
term is first rewritten 

' 1 V (B.6) 



(1 - Ac) (1 + ,tf 2 c) ft + ft Vft(l-Ac) A(l + &c)/' 

The angle can then be generated according to this distribution by generating 
the angle according to the distribution (1 — Z^c)" 1 with probability 



In 



1+01 



Pl = - T *±> _ (B.7) 



In ( j#) + In (1±|) 
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and according to the distribution (1 + /3 2 c) 1 with probability P 2 — 1 — P\. 

The full distribution can easily be generated using rejection techniques, the 
rejection weight 

( 1-/3? , 2(l+/3i/3 2 ) _ \-0j \ 

W - ^ 1 ( B - 8 ) 

(l-/3ic)(l+/3 2 c) 

is less than one. 

In practice we sometimes choose not to generate the angle according to the full 
distribution initially i.e. we postpone the latter rejection step until the event 
is generated in full. This is because the inclusion of the mass terms leads to 
a depletion of radiation in the direction of the charged particles, the "dead- 
cone" [44]. However this dead-cone can be filled by hard radiation and if the 
Pi corrections are included. In this case, if the generation of the angles is done 
according to the full distribution, the algorithm becomes very inefficient. 

In order to calculate the crude distribution we require the average photon multiplicity, 
which is given by the integral of the dipole function 



J^S( Pl ,p 2 ,k) 



= -f Z,Z 2 \n (V) ((St) ^ ( !£gjg£gj ) - 2) Ml distribution, 
= -f Z 1 Z 2 (±±ff ) In (V) ^ ( fcgjfcSj ) neglecting mass terms. 

(B.9) 

B.2 Initial-Final Dipoles 

Consider the integral of the dipole function in the rest frame of p 

J d ^S(P,P„k) = -£ Z pZl J do d*. (JL - -^)\ (B.10) 

where Z p denotes the charge on the decaying particle. As before we choose to define 
the photon momenta as being with respect to p\ hence the integral may be rewritten 

J ^ S(p, Pl ,k) = £^ Z V Z X J dc d0 dlnfc (B.ll) 

The photon energy and azimuthal angle are generated in exactly the same way as in 
appendix (|B.1| ) with the exception that now, to guarantee the possibility of conserv- 
ing momentum in the decaying particle's rest frame, the maximum allowable energy 
of the photons is 

*» = f(l-*=^'|. (B.12) 
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The generation of the polar angle 6 is more straightforward than for the final-final 
dipole. Omitting the mass terms p\ and p 2 in (|B.10|) leads to the replacement 

which may be generated by the simple mapping 



1 'l-(l + A)f^r4j I (B.14) 



The full distribution can be recovered by weighting and rejecting the events from 
this approximate distribution, with weight 

As with the final-final dipole the integral of the dipole function gives the average 
photon multiplicity for r cru d e : 

= f Z p Z, In (^) U In - 2) full distribution ( B .l 6 ) 

= ^ Z p Z 1 In (^ ai ) (-^ In ffzj§H J neglecting mass terms. 
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